## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.2.1     ✔ dplyr   1.1.2
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## 
## Attaching package: 'lubridate'
## 
## 
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
## 
## 
## corrplot 0.92 loaded
## 
## 
## Attaching package: 'ggpp'
## 
## 
## The following object is masked from 'package:ggplot2':
## 
##     annotate

Overview

This script reads in cleaned datafiles with stream temperature and invertebrate data produced by the invert_data_cleaning.Rmd and teton_data_cleaning.Rmd scripts and performs more involved analyses of stream temperature and invertebrate density/diversity over time. A few findings of interest from previous analyses are that:

* Average summer (July-Aug) stream temperatures have not changed significantly over time

* Average summer (July-Aug) air temperatures have increased over time

* Overall, invetebrate richness has not changed over time

* Invertebrate communities in snow-fed streams appear to be more variable than those in glacier-fed streams

To build on these findings, the main sections of this script are:

I. Temperature:

1. Overall trends in stream temperature over time - Time series analysis: Site-by-site and grouped by source

2. Trends in summer mean, max, min, and degree-days over time: Site-by-site and grouped by source

4. Modeling stream temperature using SNOTEL data: Site-by-site and grouped by source

II. Invertebrates:

1. Trends in richness, abundance, and density over time: Site-by-site and grouped by source

2. Modeling richness, abundance, and density using SNOTEL data: Site-by-site and grouped by source

III. Temperature and Invertebrates

1. Test for correlations between temperature metrics (mean/max/min summer; degree days): Site-by-site and grouped by source

Read in necessary files:

stream_temp <- read.csv("Temperature/cleaned_full_datasets/temps_hourly.csv") #hourly stream temp data
invert_density <- read.csv("invert_data/cleaned_csv/full_invert_densities.csv")%>%
  rename(site = Stream)#invert density data
invert_richness <- read.csv("invert_data/cleaned_csv/invert_richness.csv") #invert richness data
snotel <- read.csv("teton_snotel.csv") #snotel data
sources <- read.csv("source_info.csv")%>%
  rename(site = stream) #source info, rename for merge

### adding source info to stream temps and invert datasets: 

stream_source <- merge(stream_temp, sources, all = T)%>%
  mutate(date1 = ymd(date1)) #Add source info to hourly temperature data; re-code date as r date format
stream_source_daily <- stream_source%>%
  group_by(site, full_name, stream_code, date1, year)%>%
  summarise(t_xbar = mean(temp_c, na.rm = T), 
            t_max = max(temp_c, na.rm = T), 
            t_min = min(temp_c, na.rm = T), 
            source = unique(source)) #calculate daily mean, min, and max temps for each site
## Warning: There were 10690 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `t_max = max(temp_c, na.rm = T)`.
## ℹ In group 41: `site = "cloudveil"`, `full_name = "Cloudveil"`, `stream_code =
##   "CV"`, `date1 = 2019-09-11`, `year = NA`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 10689 remaining warnings.
## `summarise()` has grouped output by 'site', 'full_name', 'stream_code',
## 'date1'. You can override using the `.groups` argument.
density_source <- merge(invert_density, sources, all = T)
richness_source <- merge(invert_richness, sources, all = T)

I. Temperature:

2. Modeling stream temperature using SNOTEL data

First, calculate April 1 SWE and summer air temp from snotel data; add to mean temperature dataset:

#Calc means across both sites:
snotel_means <- snotel %>%
  mutate(date1 = as.Date(date, format = "%m/%d/%y"), #recode date column as R date
         airtemp_c = as.numeric(airtemp_c))%>% #recode airtemp as numeric
  group_by(date1)%>% #group by date
  summarise(swe_xbar = mean(swe_cm), #calculate mean SWE and airtemp across both stations
            airtemp_xbar = mean(airtemp_c, na.rm = T))%>%
  mutate(mo = month(date1), yr = year(date1)) #Add month and year cols
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `airtemp_c = as.numeric(airtemp_c)`.
## Caused by warning:
## ! NAs introduced by coercion
#Calculate max SWE and mean summer airtemp for each year:
snotel_summary <- snotel_means %>% 
  group_by(yr)%>% #group by year
  mutate(swe_max = max(swe_xbar))%>% #extract max SWE and add to new col
  ungroup()%>% #ungroup to get all columns back
  filter(mo == 7 | mo == 8)%>% #extract summer months
  group_by(yr)%>% #group by year
  summarise(airtemp_summer = mean(airtemp_xbar), #calculate mean summer air temp
            swe_max = unique(swe_max))%>% #keep max SWE
  rename(year = yr)

#Merge with long-term temperature averages: 
stream_snotel <- merge(summer_3plus, snotel_summary)

A. Site-by-site

i. Correlations between SWE and Summer Temp:

LMs - first, nest by site:

stream_snotel_nested <- stream_snotel %>%
  nest(data = - site)

lm for each site:

swe.temp.lms <- stream_snotel_nested %>%
  mutate(model = map(data, ~lm(summer_xbar~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "swe_max")  #remove intercept term to just get info on year for each site
swe.temp.lms
## # A tibble: 11 × 7
##    site       data              term     estimate std.error statistic p.value
##    <chr>      <list>            <chr>       <dbl>     <dbl>     <dbl>   <dbl>
##  1 mid_teton  <tibble [6 × 12]> swe_max  0.000784   0.00325     0.241  0.821 
##  2 windcave   <tibble [3 × 12]> swe_max -0.00146    0.0127     -0.115  0.927 
##  3 n_teton    <tibble [6 × 12]> swe_max -0.0561     0.0405     -1.39   0.238 
##  4 s_cascade  <tibble [3 × 12]> swe_max -0.0165     0.00937    -1.77   0.328 
##  5 s_ak_basin <tibble [4 × 12]> swe_max  0.00414    0.00603     0.687  0.563 
##  6 s_teton    <tibble [5 × 12]> swe_max -0.0487     0.0449     -1.08   0.357 
##  7 delta      <tibble [4 × 12]> swe_max -0.00178    0.00155    -1.15   0.370 
##  8 paintbrush <tibble [5 × 12]> swe_max -0.0634     0.0543     -1.17   0.327 
##  9 grizzly    <tibble [4 × 12]> swe_max -0.0567     0.0840     -0.675  0.569 
## 10 skillet    <tibble [3 × 12]> swe_max -0.0305     0.00316    -9.67   0.0656
## 11 cloudveil  <tibble [3 × 12]> swe_max -0.0166     0.0329     -0.504  0.703

Add P values and merge:

swe_summer_ps <- merge(swe.temp.lms, stream_snotel)

swe.summer.temp <- ggplot(swe_summer_ps, aes(x = swe_max, y = summer_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Annual maxium SWE, cm", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp
## `geom_smooth()` using formula 'y ~ x'

ii. Summer stream temp ~ summer air temp:

lm for each site:

air.temp.lms <- stream_snotel_nested %>%
  mutate(model = map(data, ~lm(summer_xbar~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "airtemp_summer")  #remove intercept term to just get info on year for each site
air.temp.lms
## # A tibble: 11 × 7
##    site       data              term        estimate std.error statistic p.value
##    <chr>      <list>            <chr>          <dbl>     <dbl>     <dbl>   <dbl>
##  1 mid_teton  <tibble [6 × 12]> airtemp_su…   0.0643    0.0829   0.775     0.481
##  2 windcave   <tibble [3 × 12]> airtemp_su…   0.112     0.273    0.410     0.752
##  3 n_teton    <tibble [6 × 12]> airtemp_su…  -0.457     1.32    -0.346     0.747
##  4 s_cascade  <tibble [3 × 12]> airtemp_su…  -0.251     0.0964  -2.60      0.233
##  5 s_ak_basin <tibble [4 × 12]> airtemp_su…   0.348     0.250    1.39      0.299
##  6 s_teton    <tibble [5 × 12]> airtemp_su…   0.0288    4.21     0.00683   0.995
##  7 delta      <tibble [4 × 12]> airtemp_su…  -0.214     0.0826  -2.59      0.122
##  8 paintbrush <tibble [5 × 12]> airtemp_su…  -0.0494    5.21    -0.00947   0.993
##  9 grizzly    <tibble [4 × 12]> airtemp_su…   3.18      5.68     0.559     0.632
## 10 skillet    <tibble [3 × 12]> airtemp_su…  -0.346     2.58    -0.134     0.915
## 11 cloudveil  <tibble [3 × 12]> airtemp_su…   1.81      1.04     1.74      0.332

Add P values and plot:

air_summer_ps <- merge(air.temp.lms, stream_snotel)

swe.summer.temp <- ggplot(air_summer_ps, aes(x = airtemp_summer, y = summer_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Mean daily July-Aug. air temp, C", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp
## `geom_smooth()` using formula 'y ~ x'

B. Grouped by source

First, add SNOTEL summary info to mean summer temps grouped by source:

source_snotel <- merge(summer_source, snotel_summary)
i. LMs of mean/min/max temp ~ SWE:

Re-organize data for analysis and plotting:

source_swe <- source_snotel %>%
  select(source, swe_max, t_xbar, t_min, t_max)%>% #select temperature, source, and swe cols
  pivot_longer(!c(source, swe_max), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis

Nest data by source and measure:

sourceSwe_nested <- source_swe %>%
  nest(data = -c(source, measure))

LMs for each source and measure:

source.swe.lms <- sourceSwe_nested %>%
  mutate(model = map(data, ~lm(temp~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "swe_max")  #remove intercept term to just get info on year for each site
source.swe.lms
## # A tibble: 9 × 8
##   source   measure data             term    estimate std.error statistic p.value
##   <chr>    <chr>   <list>           <chr>      <dbl>     <dbl>     <dbl>   <dbl>
## 1 glacier  t_xbar  <tibble [6 × 2]> swe_max -0.00523   0.00895    -0.585  0.590 
## 2 glacier  t_min   <tibble [6 × 2]> swe_max -0.00349   0.00961    -0.363  0.735 
## 3 glacier  t_max   <tibble [6 × 2]> swe_max -0.00689   0.00801    -0.861  0.438 
## 4 snowmelt t_xbar  <tibble [6 × 2]> swe_max -0.0491    0.0349     -1.41   0.231 
## 5 snowmelt t_min   <tibble [6 × 2]> swe_max -0.0415    0.0249     -1.67   0.170 
## 6 snowmelt t_max   <tibble [6 × 2]> swe_max -0.0548    0.0502     -1.09   0.336 
## 7 sub_ice  t_xbar  <tibble [6 × 2]> swe_max  0.0192    0.00656     2.92   0.0431
## 8 sub_ice  t_min   <tibble [6 × 2]> swe_max  0.0231    0.00566     4.07   0.0152
## 9 sub_ice  t_max   <tibble [6 × 2]> swe_max  0.0106    0.0139      0.765  0.487

Add p-values and plot:

sourceSwe_pval<-merge(source_swe, source.swe.lms)%>%
  mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position

source.swe<- ggplot(sourceSwe_pval, aes(x = swe_max, y = temp, color = measure))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source, scales = "free_x")+
   scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
  labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))+
  geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.swe
## `geom_smooth()` using formula 'y ~ x'

Save:

ggsave("Temperature/plots/source_swe.pdf", source.swe, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
ii. LMs of mean/min/max temp ~ Air temp:

Re-organize data for analysis and plotting:

source_air <- source_snotel %>%
  select(source, airtemp_summer, t_xbar, t_min, t_max)%>% #select temperature, source, and swe cols
  pivot_longer(!c(source, airtemp_summer), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis

Nest data by source and measure:

sourceAir_nested <- source_air %>%
  nest(data = -c(source, measure))

LMs for each source and measure:

source.air.lms <- sourceAir_nested %>%
  mutate(model = map(data, ~lm(temp~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "airtemp_summer")  #remove intercept term to just get info on year for each site
source.air.lms
## # A tibble: 9 × 8
##   source   measure data             term    estimate std.error statistic p.value
##   <chr>    <chr>   <list>           <chr>      <dbl>     <dbl>     <dbl>   <dbl>
## 1 glacier  t_xbar  <tibble [6 × 2]> airtem…    0.249     0.221     1.13    0.322
## 2 glacier  t_min   <tibble [6 × 2]> airtem…    0.300     0.219     1.37    0.243
## 3 glacier  t_max   <tibble [6 × 2]> airtem…    0.124     0.229     0.542   0.616
## 4 snowmelt t_xbar  <tibble [6 × 2]> airtem…   -0.701     1.10     -0.635   0.560
## 5 snowmelt t_min   <tibble [6 × 2]> airtem…   -0.712     0.806    -0.884   0.427
## 6 snowmelt t_max   <tibble [6 × 2]> airtem…   -0.488     1.53     -0.318   0.766
## 7 sub_ice  t_xbar  <tibble [6 × 2]> airtem…    0.147     0.315     0.468   0.664
## 8 sub_ice  t_min   <tibble [6 × 2]> airtem…    0.312     0.321     0.972   0.386
## 9 sub_ice  t_max   <tibble [6 × 2]> airtem…   -0.117     0.410    -0.284   0.790

Add p-values and plot:

sourceAir_pval<-merge(source_air, source.air.lms)%>%
  mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position

source.air<- ggplot(sourceAir_pval, aes(x = airtemp_summer, y = temp, color = measure))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source, scales = "free_x")+
   scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
  labs(x = "Mean July-Aug. air temp., C", y = "Stream temperature, C", color = "Measurement Type")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))+
  geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.air
## `geom_smooth()` using formula 'y ~ x'

Save:

ggsave("Temperature/plots/source_air.pdf", source.air, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'

II. Invertebrates:

2. Modeling richness, abundance, and density using SNOTEL data

A. Site-by-site

B. Grouped by source

III. Temperature and Invertebrates

1. Test for correlations between temperature metrics (mean/max/min summer; degree days)

A. Site-by-site

B. Grouped by source